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ABSTRACT 


This thesis covers the design, simulation and analysis of a SIMULINK system 
designed to predict the maneuvering characteristics of the Total Ship System 
Engineering (TSSE) program's first proposed hull design. The system is developed 
in three degrees of freedom. The ship's hydrodynamic derivatives are predicted in 
MATLAB code, while the engine is modeled completely in a SIMULINK 
environment. 

To test the system's applicability, an underway replenishment scenario is used 
to simultaneously test the steering and engine control subsystems. 

Two controllers are employed in the system. The first is used to drive the 
ship in a fashion similar to that of a human conning officer during an underway 
replenishment. The other is a root locus design used to improve the engine's 
response. 
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I. INTRODUCTION 


A. BACKGROUND 

The prediction of a ship's maneuverability as a part of the design process is becoming 
a reality. Valuable feedback can be given to a design team during the design process to 
provide insight regarding a proposed design's performance. Previously, the primary method 
of testing a hull design involved a scale model subjected to various tests in a wave tank. 
Construction and testing of models can be expensive and involve many man hours. It is 
important that once a hull's design reaches this phase that the design be as complete as 
possible. If a computer model is available for simulating a proposed hull's performance, the 
design team's productivity can be increased. The result will be a more accurate scale model 
for tank testing and, ultimately, a better ship design at lower cost. 

Toward this end, many studies have been conducted in various formats, with varying 
degrees of complexity. A low order, easy to reproduce, simplified model, coupled with a 
scenario to test a hull's design, would prove useful to design teams. This study is aimed at 
producing a simplified model of ship and power plant dynamics to be used in simulation 
studies of ship maneuverability and station-keeping. This simulation model will then be used 
to study a station-keeping control algorithm in the presence of adverse disturbances due to 
wind and wave conditions. Further, this study is aimed at providing a model for the Total 
Ship System Engineering (TSSE) program at the Naval Postgraduate School. The TSSE 
program produces, as a capstone design project, one proposed ship each year by the students 
in the program. The ships designed by the TSSE program currently have no maneuverability 
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prediction program. It is hoped that the final product of this study will provide the TSSE 
design teams with a tool enabling further study of their proposed ships form the standpoint 
of maneuverability and station-keeping. 

B. OBJECTIVES 

The goal of this study is to construct a ship-steering simulation model, complete with 
power plant dynamics and wind and wave disturbances, which will be used to assess 
alternative control laws for autonomous station-keeping during at-sea replenishment. This 
simulation will ultimately be provided to the TSSE program for use by future ship design 
teams. Specific objectives include: 

1) . The development of a mathematical model for ship dynamics, power plant and 
rudder dynamics, and disturbance factors for wind and wave conditions in various sea states. 
The relevant mathematical equations appear in Chapter II below. 

2) . The development of a computer simulation based on the equations derived in the 
mathematical model. This simulation development is described in Chapter HI below. 

3) . The development of a controller design, as described in Chapter IV below. 

4) . Test and evaluation of one controller design using the simulation. These results 
are discussed in Chapter V below. 

Input for the simulation study will be taken from the final design report for the first 
ship designed by the TSSE program: RDS-2010. 
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H. MATH MODEL DEVELOPMENT 


A. COORDINATE SYSTEMS 

Two coordinate systems are used in the description of the ship's movement relative 
to the earth. The first is a system fixed with respect to the earth oriented at some fixed point 
on the surface; the other is a system fixed with respect to the three major axes of the moving 
ship. A ship-fixed coordinate system has the advantage of constant moments of inertia and 
constant moment arm lengths in all three directions. An earth-fixed coordinate system, if 
used, would constantly change the lengths of moment arms in each direction as the ship turns 
and rolls. Convention applies control and disturbance forces to the ship in ship-fixed 
coordinates, resulting in motion of the ship’s coordinate origin relative to the earth. The 
motion of the coordinate axis of the ship represents the ship's movement with respect to the 
earth. Figure 1 illustrates the coordinate system being used. It is a right hand system 
positive forward, starboard, and down, using x, y, and z, respectively, to denote the ship 
fixed axes. Also shown are the Xo and Yo axes, which are considered fixed in the earth. 
Xo(t) and Yo(t) display a distance the ship has traveled with respect to the earth origin point 
at a time (t). The rudder angle 8 is shown as convention employs; positive rudder is for a port 
turn. The angle *¥ represents the angle the ship's velocity vector is displaced from a parallel 
of the Xo axis placed at the ship's center of gravity. For simplicity, the ship's velocity vector 
will be considered to be along its longitudinal axis, making the angle representative of the 
ship's heading with respect to the earth fixed origin (X^Y^). 
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Figure 1. Earth Centered and Ship Centered Coordinate Systems 

B. EQUATIONS OF MOTION 


An abbreviated presentation is made of a ship's equations of motion. A more in depth 
discussion of these equations can be found in Principles of Naval Architecture . 

[Ref. 2:Ch. 8] 

1. Earth Fixed Equations 

From Principles of Naval Architecture [Ref 2], a ship's movement with respect 
to an earth fixed reference point (X^Y^) can be described by equations (la and b) as: 

(a) X 0G = ucosy - vsiny 

( b ) Y og = vcos\j/ + wsinvj/ 
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where: x, y, and z are in ship fixed coordinates, and u = x, v = y, \j/ = ship's heading, 
are termed surge, sway and yaw respectively. 

2. Ship Fixed Equations 

If a moving set of axes is desired, Ref 2 shows that a ship's equations of 
motion on a horizontal plane are: 

(P) X u (w - u x ) + (m - X.) u = 0 

(b) -Y v v + ( m -Y.)v- (Y r - m) r - (Y r - mx G ) f = Y s 5 (2) 

(c) -Nv - (N. - mx G ) v - (N r - mx Q ) r + (J 2 - N.) f = N s 5 

where: X = total force in the x direction 
Y = total force in the y direction 
N = total moment about the z axis 
I 2 = mass moment of inertia about the z axis 
m = ship's mass 

Xg = distance ship's center of gravity is displaced 
from the centerline (positive forward) 

8 = rudder deflection angle 

r = 'F 

Equations (2a-c) assume symmetry about the ship's longitudinal axis, and employs 
a shorthand notation where: 


5 




(a) Y v = dY/dv 

(b) N. = dN/dr, etc.... 


( 3 ) 


Assuming linearity of the above, we see that X in equation (2a) is not coupled with 
second two equations (2a and b), hence it will be dealt with separately. 

a. Nondimensionalizing Variables 

Equations of (2a and b) has units offeree, while equation (2c) has units of a 
moment. Therefore, it is convenient to nondimensionalize the variables as follows: 


(«) 

0 ) 

(c) 

(d) 

(e) 

if) 

ig) 
0 ) 
(0 
if) 


_ / _ 
TYl = 


m 


p/2 V 
jf m 

2 p/2 L 5 

\7 

yL = 


p/2 L 7 U 2 


Y' = 


Yi = 


y! = 


;v' = 


JV/ = 


K = 


n! = 


p/2 Z 3 f/ 

n 

p/2 VU 

h 

p/2 L A U 

K 

p/2 L 3 U 2 

K 

p/2 Z 4 Z7 

p/2 X 4 t/ 
N 


p/2 L 


(4) 
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where: p = density of seawater 


L = ship's length 
U = ship operating speed 

b. Formulae for Predicting the Hydrodynamic Derivatives 
From Regional Deterrence Ship (RDS-2010) [Ref. l:Ch. 8: Sec. 10], the 
following formulae are found for deriving hydrodynamic derivatives from hull data: 

(«) K -- (rX * (Xl), 

(b) r; .<r\ 1/2( 7 v V 

W ri - (XX + 04o4 

(d) y. = 0 - 1 / 2(74 

(«) ti-wX i mxi>, 

W if! = <K\ ♦ i/4(r'y (S) 

(g) ^ = Qf\ * 1/4 VIdot), 

( h ) n ! = o -1/2(74 

(0 

CO A/i, = -1/2(74 

where: () h denotes the contribution to the variable from the hull, and () f denotes the sum 
of contributions to the variable from fins and other appendages as shown in equations (6a-g). 
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where: 


id) (Y\ = -UTIL + C D 
(b) (YX = k x m[ + xJL{Yl) h 
ip) QfX = -Q»t -kjjn) + xJUXX 
id) QtX = m' z x/L + (#4 

5^/72 

(*) (i’/)* = f («) 

bow 

stem 

W ^ = * 7 7 ^ / c / 2 * 2 ^ 

^ * bow 


(*) 

(*) 

(*) 


stern 



m' z = Qc'/k^) m 2 


stem 

f C s d 2 xdx 

bow 

stem 

/ c ^ Ia 

dew 


(7) 
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and as empirically derived by Vann [Ref. 3 ]: 


(a) k x = 0.3(2 TIL) 

(b) k 2 = 1.0 - 0.5(2 TIL) 

(c) k' = 1 - 1.333(2 TIL) 

id) C s = C Q + C x (sa/bd) + C 2 (salbd) 2 (8) 

(e) C 0 = 0.8572 + 0.5330(4^/6) 

(f) C, = 3.374 - 1.3661(4^6) 

(g) C 2 = -1.7323 + 0.8670(4^ 

where: 

T= ship's draft 

C D = ship's coefficient of drag 
kj = longitudinal coefficient of accession to inertia 
k 2 = lateral coefficient of accession to inertia 
k’ = rotational coefficient of accession to inertia 
C s = two dimensional sectional inertia coefficient, calculated using strip 
integration along the ship’s hull, 
sa = section area of section being considered 
d = draft at section being considered 
b = beam at section being considered 
Co, C„ C 2 = interim variables used in the calculation of C s 
Xp = distance from centerline to point of application of fluid force 
x^L = Vi prismatic coefficient 
Af = profile area of the appendage 
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Af' = nondimensionalized profile area of fin or appendage (A/LT) 
a = aspect ratio of the fin or appendage, defined by Gillmer [Ref. 4] as: 


a = 



( 9 ) 


h = span of the appendage 

Calculation of the above derivatives is accomplished in MATLAB code with the files 
HYDROGEN.M and CSGEN.M that are listed in appendices A and B. 

3. Propulsion Equations 

To solve for motion along the ship's longitudinal axis, equations (10a and b) for gas 
turbine engines are used, as presented by Tozzi. [Ref. 5] 

(a) V = g/m(T p - R s ) 

(i) JV = -Ufi, - Q, - Q p ) (10) 


where: V = ship velocity in feet per second 

N = Propeller shaft speed in rotations per second 

g = acceleration due to gravity in feet per second squared 

m = ship mass in pounds mass 

T p = propeller thrust in pounds force 

Rs= resistance of the ship in pounds force 

I = moment of inertia of the drive train referred to the propeller shaft 
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Q t = engine torque available at the propeller shaft in foot pounds force 
Q f = friction torque losses in foot pounds force 
Q p = propeller torque losses in foot pounds force 

a. Propeller Thrust 

The thrust is a measure of force developed by the propeller in the direction of 
the shaft. Mathematically, it is expressed as: 

T P = k T p D 4 N 2 (11) 


where: k x = thrust coefficient 

p = density of seawater in slugs per cubic foot 
D = propeller diameter in feet 


b. Ship Resistance 

The resistance a hull form has is normally found during scale model tests in a 
wave tank. For the RDS -2010 the resistance is estimated as a function of speed. From data 
provided by Alexander [Ref. 1] it is seen that ship resistance varies as: 


R s = 172 V 2 


( 12 ) 


c. Drive Train Moment of Inertia 

As presented by Tozzi [Ref. 5], the moment of inertia includes the 
contributions of all drive train components referred to the propeller shaft. For this simulation. 
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Equation (13) is used to estimate the moment of inertia, based on Ref. 5. 


/ = 2.0 x 10 5 


( 13 ) 


d. Shaft Torque 

The shaft torque (Q x ) as supplied in the report on the RDS-2010 includes the 
effect of the reduction gears in the drive train, and is a measure of the torque, in foot pounds 
force, provided to the propeller shaft. 

e. Friction Torque 

Friction in the entire drive train can be modeled as a loss in torque developed 
to produce thrust. From Ref 2 it is seen that friction torque can be approximated as: 

Q f = 6000 N (14) 


f. Propeller Torque 

The torque required to rotate a propeller under various conditions is normally 
measured during open water tests. For the subject ship it can be calculated using data 
provided in Ref. 1 as: 

Q p = K q p N 2 D s (15) 

where: Kq = Torque coefficient for the hull under consideration. 
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4. Steering Equations 

A low order rudder model can be developed similar to Ref. 6. It is modeled as an 
integrator with gain Kg' as follows: 


(«) 

(*) 




u 


K = 


(16) 


where: 5 m = maximum rudder angle 
8einax = maximum error input. 


5. Ship Dynamics 

To convert the Y and N equations of Equations (2a and b) to a form usable in the 
simulation it is necessary to develop matrices as follows: 


(a) [M] x = [H] x + [R] u 

(b) x = [M\~ l [H] x + [AC 1 [^] « 


(17) 
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where: 


(«) 

(b) 

(c) 

id) 

00 


m 

m 

m 


im' - Yfom'x' - y!) 
(w'x' - <)(// - N!) 

[ Y'(Y; - m 0 

[iVvW 


x = [v r] r 
u = 3 


( 18 ) 


The new matrices formed are implemented as follows: 


id) 

x = Ax + Bu 


ib) 

g 

**H 

5 

u_Li 

11 

(19) 

(c) 

b = m -■ [A] 



C. DISTURBANCES 

To test the system's performance under varying conditions, wind and sea disturbances 
are added to model the effect of increasing sea state on the system. Also included is the 
venturi effect, a phenomenon affecting the sway forces and yaw moments of ships as they pass 
each other in close quarters. 

Lastly, measurement noise is added to the observation of the supply ship's position to 
model the effect of winds and seas acting on the supply ship. This noise models the effect of 
sea state by increasing position estimate error as sea state increases. 
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1. Sea Disturbances 


In order to simply model what can be made infinitely complex, Uhrin [Ref. 7] employs 
formulas (20a-f) for a wave. These formulas generate a disturbance at a primary frequency 
and its second harmonic: 

(a) W = WF (1 + WRV) sin (WEF) + s i n (2 WEF) 

(b) mF = ink (u + ws cos (XPrs) 

(c) - V SHIp (20) 

(d) X s - W costae 

(e) Y s = W sinT^ 

if) N s = W sin (2 

where: W = total wave force 

WEF = wave encounter frequency 

WF = maximum wave force 

WL = wave length, normalized by WLT/LOA 

WLT = true wave length 

LOA = ship length overall 

WS = wave speed in feet per second 

Trs = relative seas heading 

T ts = true seas heading 

'J'shp = ship heading 

X,, = sea force in x direction 
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Y s = sea force in y direction 
N s = sea moment in N direction 
WRV = wave random variable 

The wave random variable is added to model the randomness of the sea. It is 
implemented as zero mean gaussian white noise with a variance of 0.01. 


2. Wind Disturbances 

To account for the wind's effect on the ship, a model is designed similar to that 
employed by Clark. [Ref. 8] It is necessary to take wind true velocity and speed in the earth 
fixed reference frame, and translate each component of the forces and moment to the ship 
fixed reference frame. 

First, the relative velocity and direction must be found: 

( a ) Vrw 

0) \f RJfr 

V w COS VsHIp) u 


Ml ~ 


( 21 ) 


Coefficients for each component of wind effect are developed as follows: 


(«) 

0 ) 

(c) 



Pa A f 

7000p w LBP 2 
- P ajs 
8000p w LBP 2 

Pa £s 

p. LBP* 


( 22 ) 
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Lastly, the nondimensional forces and moment are developed in equations (23a-c): 

(a) X!r = (-^) 2 C m sin[9/7(| Vw | - n)] 

(*) sm(y, RW ) 

(<0 < = sin(2^) 

where: X w ' = wind force in x direction 
Y w ’ = wind force in y direction 
N w ' = wind moment about yaw axis 
Yrw = relative wind velocity 
^rw = relative wind direction 
V w = true wind velocity 
V P TW = true wind direction 
'f'sHjp = ship heading 
u = ship surge 
v = ship sway 

Cwx = x direction wind coefficient 
Cwy = y direction wind coefficient 
Cwn = N direction wind moment 
p a = mass density of air 
p w = mass density of water 


( 23 ) 
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Af = ship frontal area 

A* = ship side area (estimated as 10Af) 

LO A = ship length overall 

LBP = ship length between perpendiculars 

X w ' = nondimensional x direction wind force component 

Y w ' = nondimensional y direction wind force component 

N w ‘ = nondimensional N direction wind force component. 

3. Venturi Disturbances 

The venturi effect can be closely estimated as a nonlinear function of 
longitudinal separation that is inversely proportional to lateral separation as seen in Ref. 6. 
Equations (24a and b) are used to approximate this: 

(a) Y'y^ = 0.45) X 10- 5 [smc (0.0029 IS)] Q 

„/ _ _ „/ < 24 ) 

\°) iy VENT ~ 1 VENT 


where: 

Yvent' = nondimensional venturi force 
Nvent' = nondimensional venturi moment 
LS = longitudinal separation in feet 

Q = lateral separation multiplier (varies linearly from 1 to 0.5 as separation increases 
from 100 to 200 feet). 
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4. Measurement Noise 


The measurement noise, as described above, is simply modeled as gaussian 
white noise with a variance of 5, 10 or 20, which will vary with sea states 0, 3, and 6, 
respectively. 
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m. SIMULATION DEVELOPMENT 


A. PROBLEM GEOMETRY 

The underway replenishment scenario considered in this research will place the supply 
ship at an earth fixed origin (Xo,Yo) at problem start. Throughout the problem the supply 
ship will remain at the relative origin. The receive (subject) ship will start the problem in 
waiting station, one thousand yards astern of the supply ship. Both ships will begin the 
problem on base course, at base speed (course 000 at 14 knots). This system will drive to a 
starboard/ port supply/ receive geometry. The system will seek to drive the receive ship to 
a position one hundred and thirty feet abreast of the supply ship, at (0,130) relative position. 



Ysd =y separation desired =130 

_ ......._i 


I 

1-1 Yr(t) = 130-yr(t): y remaining 

S(t) 

4(0 


\ y(t) / D (0 ~ sqrt[xrft) A 2 + (130-yr(0) A 2] 

So 

/ Ao 

Do = sqrt[3000 A 2 + 130 A 2] 


Vy>^R(0 = (xrft),yr(0) 

Ro = (-3000,0) 


Figure 2. Problem geometry. 


The initial positions of the supply and receive ships is displayed in Figure 2 as S 0 and 
Ro respectively. The initial aim point is denoted as A„. From the figure it is evident that the 
aim point moves with the supply ship at base course and speed. As the problem progresses. 
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the ship is at a relative position (xr(t),yr(t)), steers an angle v|/(t), and has a distance to travel 
D(t). Thus, the initial angle vj/ 0 and distance D 0 are not constantly fixed with respect to the 
absolute origin, and must be computed and fed into the system as command signals for all 
time (t). 

B. SHIP'S PLANTS 

1. Rudder 

As previously mentioned, a low order model for a rudder plant was designed in Ref. 
6, and is shown as Figure 3. Heading error is input to the system at the in block. The rudder 
stops limit the amount of travel and is modeled as a saturation block. The rate limit block 
simulates a hydraulic pump that limits the amount of error into the rudder, which is modeled 
as an integrator with gain Kg’. The rudder angle (6) is the output of the subsystem at the out 
block. The SIMULINK implementation requires renaming the variable Kg' in the rudder 
subsystem as Kgpri. 



Figure 3. Rudder subsystem. 
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2. Ship's Engine 

The Engine model is designed to output the thrust necessary to drive the ship. Input 
is the desired torque. Next, the subsystem develops the shaft acceleration as defined earlier. 
An integrator is used to generate the propeller rotations per second, which is then converted 
to a thrust value. The initial condition set on the integrator starts the problem with the ship 
already at base speed. The thrust output is then calculated as a function of shaft rotations 
as noted previously. From Equations (10a and b) the conversion into a SIMULINK 
environment is relatively straightforward as seen in Figure 4. 



Figure 4. Engine Subsystem. 
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C. SHIP’S DYNAMICS 


1. Steering Dynamics 

To model the ship dynamics in SIMULINK, the matrices in Equations (19a-c) are 
used in the subsystem. These matrices lend themselves very well to a SIMULINK 
implementation. The subsystem also adds the Y and N disturbances to the accelerations. 
Figure 5 displays the ship dynamics subsystem. 



Figure 5. Steering dynamics subsystem. 


2. Propulsion Dynamics 


Equation (10A) is used to produce the ship's longitudinal velocity from the developed 


propeller thrust. Again, the conversion is uncomplicated as seen in Figure 6. 
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Figure 6. Propulsion dynamics subsystem. 


D. EARTH ORIENTED POSITION GENERATION 


Equations (la and b) are very simply translated into a SIMULINK subsystem to solve 
for velocities in the Xo and Yo directions. These velocities are integrated to obtain the 
positions with respect to the earth fixed origin. Figure 7 shows the earth position generation 
subsystem as implemented. 
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Figure 7. Earth oriented position generation subsystem. 


E. RELATIVE POSITION GENERATION 

In order to generate the ship's position relative to the supply ship, the position of the 
supply ship with respect to the earth fixed origin must be compared with the receive ship's 
position with respect to the same origin at each time interval. This is accomplished by using 
base speed as the ideal supply ship motion along the Xo axis. The supply ship's motion along 
the Yo axis is ideally zero. To both position estimates of the supply ship, the measurement 
noise is added. This effectively models a ship's sensor (radar), and accounts for the supply 
ship's motion in an increasing sea state. The relative position generator subsystem is shown 
in Figure 8 as implemented. 
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Figure 8. Relative position generation subsystem. 


F. DISTURBANCE GENERATION 


As discussed in the previous chapter the disturbances are to test the system's 
performance under less than ideal conditions. The forces and moments generated are summed 
from each component. 
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1. Wind Disturbance 


From the previous chapter the Figure 9 shows the subsystem developed to generate 
wind disturbances. 



Figure 9. Wind disturbance subsystem. 
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2. Sea Disturbance 


The mathematical model for seas as developed previously translates very well into a 
SIMULINK subsystem, and is shown in Figure 10. 
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3. Venturi Disturbance 


The venturi phenomenon is modeled in SIMULINK with the subsystem shown as 


Figure 11. The subsystem follows directly from Equations (27a and b). 



Figure 11. Venturi disturbance subsystem. 
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IV. CONTROLLER DESIGNS 


Two types of controllers are employed within the system designed. The first is a 
command guidance signal styled to drive the ship similar to a conning officer during an 
underway replenishment. This is accomplished by the conversion of the receive ship's position 
relative to the supply ship into command signals to the rudder and engine. The second is a 
root locus design for the engine plant to drive the ship more efficiently. 

A. COMMAND SIGNAL GENERATION 

This first controller design is implemented as a separate subsystem, converting the 
receive ship's relative X and Y positions into rudder and engine command signals as 
mentioned above. In order to achieve the desired goal of imitating a conning officer, the first 
goal is to have the ship open the lateral separation in order to safely close the longitudinal 
distance without danger of collision. While the rudder steers the ship towards the aim point, 
the engine must be sped up to close the longitudinal separation. A typical conning officer will 
bring the ship up to a speed of 20 knots to expeditiously gain station. Once the ship is within 
approximately 200 feet of station, the ship is slowed to 1 knot above base speed to avoid 
going past the station. 

Once the ship is in station, the goal is to stay there. This requires small, continual 
adjustments in both the rudder and the engine in response to the effects of the wind, sea, and 
the presence of the supply ship. 

To gain and then maintain station defines two separate problems as defined above. 
The first is an approach problem requiring large corrections. The second is a station-keeping 
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problem requiring finer adjustments to maintain the desired position and avoid dangerous 
maneuvers while in close quarters to the supply ship. 

During the approach phase, the desired heading is conveyed to the rudder subsystem 
as 'P(t), as defined previously. The distance D(t) is used to generate a velocity command 
signal ordering the engine to 20 knots until within 200 feet of station. The velocity command 
signal is also converted to a torque command signal based on the predicted torques at various 
speeds for the subject ship. 

The lateral separation, being opened to the proper distance, is employed as the 
division between approach and station-keeping phases. Once the lateral distance is within 5 
feet of the desired setting, the command generation subsystem toggles to the station-keeping 
phase. During this phase the rudder is restricted to a maximum amplitude of one-fifth of 'F 0 . 
The rudder will steer the ship with minute adjustments when within 20 feet of the desired 
lateral separation. The station keeping phase also restricts the engine response by basing the 
D(t) on only the longitudinal separation. Further, the velocity command is regulated to within 
1 knot of base speed when within 40 feet of station. The engine will respond to gradually 
speed up or slow down in response to the relative longitudinal position. The rather 
complicated verbal description above is actually very simple to implement in SIMULINK. 
The relative positions are used as input and converted to 'f'(t) and D(t) signals as described 
in the problem geometry. These are used with look-up tables to accomplish the intent of the 
subsystem. The subsystem is pictured in Figure 12 as implemented. 
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Figure 12. Command signal generation subsystem. 


B. ROOT LOCUS CONTROLLER 

The engine requires a separate controller developed to improve the response of the plant. 
Figure 13 displays the actual plant implemented. 
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Y = u 



Figure 13. Actual engine plant. 
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where: 


<«) 


S X = 


u - u Q 
N-N 0 


(P) 5 Q t ^Q t - Q Tq 


(27) 


The first term in the above non-linear Equation (27b) is assumed to go to zero, leaving 
a linear equation. The second term defines the plant's A matrix. The third term defines the 
plant's B matrix. For the controller design, a new 5Q T is defined as below. 

5 Q t = K (u - w 0 ) = K Sv (28) 


Then for the approximate system: 

5x = A 8x + B K 5u = A fix + B K C fix (29) 


hence: 


(a) 5x = [A - KB C\ hx 

(b) 5y = C 5x 

because K is a scalar. 

The poles of this system are at: 

DET(sI - A + KBC) 


(30) 


(31) 


and, the system transfer function is: 

C [si - A = + KBC ]' 1 B 


(32) 
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A suitable value for K was found, experimentally to be 100000. This was employed 


in the system as displayed in Figure 14. 



Desired 

Velocity 

Input 


Figure 14. Engine plant as implemented. 
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V. SIMULATION RESULTS 

The system was run several times to demonstrate its performance under varying 
conditions. First the system was run without disturbances as shown in Figures 15-17.. Next, 
the system was run with only the measurement noise applied as a disturbance seen in Figures 
18-20. For this run the measurement noise was increased to 20 to better display its effects 
on the system. Lastly, the system was run with wind, sea, and venturi disturbances applied 
in addition to the measurement noise shown in Figures 20-23. For the last run the wind and 
seas were set at 5 knots, the normalized wavelength set at 0.1, and the true direction set at 
000 degrees true (on the bow). For each run three plots are shown. The first plot presents 
the relative trajectory of the receive ship (a bird's eye view of the track taken by the receive 
ship). The next two plots display the lateral and longitudinal trajectories plotted against time 
to compare settling times of the different cases. 



Figure 15. Undisturbed system relative trajectory 
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Figure 17. Undisturbed system longitudinal trajectory. 
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Figure 19. Measurement disturbance lateral trajectory. 
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Figure 20. Measurement disturbance longitudinal trajectory. 
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Figure 22. Full disturbances lateral trajectory. 
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VL CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

From the simulation results in the previous chapter it is evident that the system 
performs as desired when undisturbed. The lateral separation is opened safely while closing 
the longitudinal distance. And, once the station is attained, the system keeps the ship at the 
desired point. Additionally, it is shown that the measurement noise has little effect on the 
system. In each of these first two cases the system has a time to station of 420 seconds. 

However, it is also evident that the system does not comfortably open the lateral 
separation while being exposed to the wind, sea, and venturi disturbances. But, the system 
still has a satisfactory time to station of 430 seconds, which is very comparable to the 
undisturbed case. And, once on station the system still keeps the ship in place. 

B. RECOMMENDATIONS FOR FURTHER STUDY 

Based on the above conclusions, the first recommendation is to design a controller 
to better handle the rudder during exposure to disturbances. Further, many different control 
schemes for both the engine and the rudder may be designed and compared using classical 
as well as optimal control approaches. 

As the TSSE program has designed two additional ships since the RDS-2010, it may 
be desired to apply this study to the other ships. This may require the development of 
different propulsion plants which could be stored and modified for future use by the TSSE 
program. Additionally, the TSSE program could employ this system to test more 
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fundamental performance measures of their proposed ships, such as tactical diameter. It 
should also be noted that this system would be applicable to a commercial vessel design. 

Finally, the disturbance model would prove useful to any design team wishing to test 
the sturdiness of a proposed hull design. The forces and moments the hull would be exposed 
to can be measured, and compared to the maximum levels the structure should be exposed 
to. 
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APPENDIX A HYDROGEN.M 


% HYDROGEN.M: FILE TO GENERATE HYDRODYNAMIC DERIVATIVES, WIND 
% COEFFICIENTS, WIND AND SEAS TRUE DIRECTION, AND ENGINE 
% SIMULATION CONSTANTS FOR SIMULINK FILE "SHIPSIM". 

% VALUES REQUIRING USER INPUT FOR SYSTEM CUSTOMIZATION ARE 
% NOTED AS (input) IN THE COMMENT LINE. 

% THIS FILE MUST BE RUN PRIOR TO "SHIPSIM" TO GENERATE THE 
% NECESSARY INPUTS 


rho=1.9914*0.0311;% Density of H20 in Slugs/Cub. Ft. 

% at 55 deg F * (Slugs/Lb Mass) 

U = 14* 1.689 ;% Nominal Ship's FWD Speed 

% (14 knots) in Feet per Sec. 

% HULL COEFFICIENT OF DRAG AT ZERO ANGLE OF ATTACK CALCULATION 

rhos= 1.9914 ;% Density of H20 

vs= 1.3535 ;% Kinematic Viscosity of H20 

q=.5*rhos*vs A 2 ;% Dynamic Pressure of H20 
Do=123559 ;% Drag Force at 16 Knots 
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WSA=24204.1 ;% Wetted Surface Area at Displacement 

Cdo=Do/(q*WSA) ;% Hull Coefficient of Drag 

% NECESSARY INPUTS FOR HYDRODYNAMIC DERIVATIVE GENERATOR 


D =5721.8 

;% (input) Ship Disp. in Long Tons (5721.8) 

T = 15.01 

;% (input) Ship's Draft in Feet (15.01) 

L =409.315 

;% (input) Ship’s Length (LOA) (409.315) 

B = 55.09 

;% (input) Ship's Beam in Feet (55.09) 

Cb =0.599 

;% (input) Ship's Block Coefficient (0.599) 

mship = (2240*D)/32.2 ;% Ship's Mass Calculator 

delmaxdot = 2 

;% (input) maximum rudder rate 

delemax = 7 

;% (input) maximum error input 

Kg = delmaxdot/delemax ;% Rudder Gain Calculator 

Xg =-5.85 

;% (input) Estimate of Xg (LCB for Disp.) 

xp =-22.5 

;% (input) Long. Center of Float. (LCF) 


% at Disp. 

Priz= 0.650 

;% (input)Prismatic Coefficient 


% APPENDAGE DATA: CALCULATED FOR RUDDERS AND THE SONAR DOME, 
% MAY BE CUSTOMIZED FOR OTHER APPENDAGES BY ADDING THE 
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NECESSARY 


% CALCULATIONS SIMILAR TO THOSE BELOW. 

ProfArudl = 117.4 ;% (input) Rudder Profile Area (117.4) 

PAnondl = ProfArudl/(L?T) ;% Non Dimensionalized Profile Area 

hrudl =11.95 ;% (input) Rudder Height 

arudl = hrudl A 2/ProfArudl ;% (input) Rud Eff Aspect Ratio 

hrudndl= hrudl/L ;% Non Dimensionalized Height 

Cdrudl = 0.008 ;% (input) Rudder’s Coef of Drag (0.008) 

mrudl = (2*pi)/(l+2/arudl);% (input) Rudder's Lift Curve Slope Calc 

ProfAmd2 = 117.4 ;% (input) Rudder Profile Area (117.4) 

PAnond2 = ProfArud2/(L*T) ;% Non Dimensionalized Profile Area 

hrud2 = 11.95 ;% (input) Rudder Height 

arud2 = hrud2 A 2/ProfArud2 ;% (input) Rud Eff Aspect Ratio 

hrudnd2= hrud2/L ;% Non Dimensionalized Height 

Cdrud2 = 0.008 ;% (input) Rudder's Coefficient of Drag (0.008) 

mrud2 = (2*pi)/(l+2/arud2);% (input) Rudder’s Lift Curve Slope Calc 

ProfAdom = 1400 ;% (input) Sonar Dome Profile Area (1400) 

PAnond = ProfAdom/(L*T) ;% Non Dimensionalized Profile Area 

adorn = 0.35 ;% (input) Sonar Dome Aspect Ratio 

hdom = 22.26 ;% (input) Sonar Dome Height 
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hdomnd= hdom/L 


;% Non Dimensionalized Height 


% CALL TO INERTIAL COEFFICIENT GENERATOR 
csgen 


% DERIVATIVE GENERATOR FOLLOWS EQN'S FROM REF. 1 CHAP-8 pp525 
% EQN’S 67 FOR HULL WITH SMALL DEADWOOD 

% NON-DIMENSIONAL HULL HYDRODYNAMIC DERIVATIVE GENERATOR 
mpri = mship/(0.5*L A 3*rho); 

Yvh = -pi*(T/L)*Cdo; 

Yrh = -(kl *mpri) + ((xp/L)*Yvh); 

Nvh =-(m2-(kl*mpri)) + ((xp/L)*Yvh); 

Nrh = -(mz*(xbar/L)) + ((.5*Priz) A 2*Yvh); 

Yvdoth = -m2; 

Yrdoth = 0; 

Nvdoth = 0; 

Nrdoth = -((kprime*pi)/(L A 3*T A 2))*int3; 

Izh = mship/(0.5*L A 5*rho); 

Kgpri = Kg*(L/U); 

Xgpri = Xg/L; 
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% NON DIMENSIONAL RUDDER (CONTROL) HYDRODYNAMIC DERIVATIVE 
GENERATOR 

Ydel =(4/3)*PAnondl *((2*pi)/(l+(2/arudl))); % 4/3 ACCOUNTS FOR TWO RUDDERS 
Ndel = -.5*(Ydel); 


% NON DIMENSIONAL APPENDAGE HYDRODYNAMIC DERIVATIVE 

Generator 

Yvfirl = -PAnond 1 *((2*pi)/( 1 +(2/arud 1))); 

Yvfr2 =-PAnond2*((2*pi)/(l+(2/arud2))); 

Yvfd = -PAnond* ((2*pi)/( 1 +(2/adom))); 

Yvap = Yvfrl+Yvff2+Yvfd; 

Yrfrl =-.5*Yvfrl; 

Yrfr2 =-.5*Yvfr2; 

Yrfd =-.5*Yvfd; 

Yrap = Yrfr 1 + Y rfr2+Yrfd; 

Nvfrl = Yrfrl; 

Nvfr2 = Yrff2; 

Nvfd = Yrfd; 

Nvap = Nvfir 1+Nvff2+Nvfd; 

Nrfrl = .25*Yvfrl; 

Nrfr2 = ,25*Yvfr2; 
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Nrfd =.25*Yvfd; 


Nrap = NrM+Nrff2+Nrfd; 

Yvdotfrl = -4*pi*PAnondl*hrudndl/sqrt(arudl+l); 
Yvdotfr2 = -4*pi*PAnond2*hrudnd2/sqrt(arud2+l); 
Yvdotfdm = -4*pi*PAnond*hdomnd/sqrt(adom+l); 
Yvdotap = Yvdotfrl+Yvdotfr2+Yvdotfdm; 

Yrdotfrl = 2*pi*PAnondl*hrudndl/sqrt(arudl+l); 
Yrdotfir2 = 2*pi*PAnond2*hrudnd2/sqrt(arud2+l); 
Yrdotfdm = 2*pi*PAnond*hdomnd/sqrt(adom+l); 
Yrdotap = Yrdotfrl+Yrdotfr2+Yrdotfdm; 

Nvdotfrl = Yrdotfrl; 

Nvdotfr2 = Yrdotff2; 

Nvdotfdm = Yrdotfdm; 

Nvdotap = Nvdotff l+Nvdotfr2+Nvdotfdm; 

Nrdotfrl = .25* Yvdotfrl; 

Nrdotff2 = ,25*Yvdotff2; 

Nrdotfdm = .25*Yvdotfdm; 

Nrdotap = Nrdotfrl+Nrdotfr2+Nrdotfdm; 

% TOTAL HYDRODYNAMIC GENERATOR 
Yvtot = Yvh + Yvap; 
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Yrtot = Yrh + Yrap; 

Nvtot = Nvh + Nvap; 

Nrtot = Nrh + Nrap; 

Yvdot = Yvdoth + Yvdotap; 
Yrdot = Yrdoth + Yrdotap; 
Nvdot = Nvdoth + Nvdotap; 
Nrdot = Nrdoth +Nrdotap; 

% MATRIX GENERATION 
al 1 = mpri-Yvdot; 
al2 = (mpri*Xgpri)-Yrdot; 
a21 = (mpri*Xgpri)-Nvdot; 
a22 = Izh-Nrdot; 

bll = Yvtot; 
bl2 = Yrtot-mpri; 
b21 = Nvtot; 

b22 = Nrtot-(mpri*Xgpri); 

ell = Ydel; 
c21 = Ndel; 
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A = [all al2;a21 a22]; 
B = [bllbl2;b21b22]; 
C = [ell ;c21 ]; 


Apri = inv(A)*B; 

Bpri=(inv(A)*C); 

Cpri = [1 0;0-l]; 

Dpri=[0;0]; 

% WIND COEFFICIENT GENERATOR FOR POOR THESIS 
% CONSTANTS 

rhow = 64 ;% Density of H20 in lbs/ft A 3 

rhoa = 0.0752 ;% Density of Air in lbs/ft A 3 

LBP = 390.00;% (input) Ship Length Between Perpendiculars 

LOA = 409.31;% (input) Ship Length Overall 

Af = 2037.9 ;% (input) Ship Frontal Area 

As = 10*Af ;% (input) Ship Side Area (Estimate) 

q=rhoa/rhow; 

% COEFFICIENT GENERATION 
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Cwx = q*Af/LBP A 2*7000 ;% Wind Surge Force Coefficient 

Cwy = q*As/LBP A 2*8000 ;% Wind Sway Force Coefficient 

Cwn = q*As*LOA/LBP A 3*1000;% Wind Yaw Moment Coefficient 


% ENGINE CONSTANTS AND TORQUE GENERATION 
1=1.9e5 ;% Prop shaft moment of inertia 

gm=(32.2 A 2/(5721 *2240)) ;% Equivalent form of g/m 

Qt=219597 ;% Input torque at 16 knots 

Kl=2*gm*0.185*1.9914*15.5 A 4 ;% Tp, Prop Thrust Calc. 

K2=gm*171.7658 ;% Rs, Ship Resistance Calc. 

K3=(6000/(2*pi*I)) ;% Qf, Fric. Torque Loss Calc. 

K4=(0.0417*1.9914*15.5 A 5/0.985)*(l/(2*pi*I));% Qp, Prop. Torque Loss Calc. 
K5=l/(2*pi*I) ;% Qt Multiplier 

Aeng=[-2*K2*27.024 2*K1*1.72;0 (-2*K4*1.72)-K3J; 

Beng=[0 K5]'; 

Ceng=[l 0]; 

Deng=[0]; 

[numeng,deneng]=ss2tf(Aeng,Beng,Ceng,Deng); 


53 



%%% torquefit 

Torq = [219597 586322 689169];% Shaft Torques at Vel, Speeds From 
Vel = [16 25.26 26.49] ;% RDS-2010 Final Report 

Speeds = linspace(16,22,60) ;% 

Torques = interpl(Vel,Torq,Speeds,'linear');% 

%plot(Speeds,T orques) 

vin = [16 17 18 19 20 21]; 

tout = interpl(Vel,Torq,vin, , linear'); 
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APPENDIX B CSGEN.M 


% FILE TO GENERATE kl,k2,kprime, AND Cs FOR RDS 2010 
kl=0.3*(2*T/L); % Empirical Formulae 

k2=l-(.5*(2*T/L)); % From Vann Thesis 

kprime= 1-(1.33*(2*T/L)); % 


% SECTION AREA AT EACH STATION 

secar=[0 0 0 25.69 52.5 158.7 271.48 382.65 485.27 574.11... 

646.09 700.07 736.36 755.95 759.89 748.75 722.26... 

668.83 585.81 470.15 324.39 164.23 30.00]; . 

% LOCAL DRAFT AT EACH STATION 

d= [0 0 0 9.611 15.006 15.006 15.006 15.006 15.006 15.006... 

15.006 15.006 15.006 15.006 15.006 15.006 15.006 14.743... 

13.861 12.187 9.516 5.716 1.308]; 

% LOCAL BEAM AT EACH STATION (PAD W/1 'S TO AVOID DIVIDE BY 0 ERROR) 
b= [0 0 0 2.452 4.82 14.418 24.208 33.062 40.386 45.984 49.934... 

52.488 53.978 54.734 55.032 55.050 54.844 54.196 52.798 50.206... 

45.996 39.996 32.574]; 
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% DISTANCE x FROM CL AT EACH STATION 
x= [204.655 194.997 185.34 179.93 174.521 154.704 134.887 115.069... 
95.252 75.435 55.618 35.801 15.984 -3.834 -23.651 -43.468... 

-63.285 -86.848 -110.410 -130.973 -157.535 -181.098 -204.66]; 

for i=l:length(d) 
if d(i) = 0 
mult(i) = 0; 
else 

mult(i) = secar(i)/(b(i) *d(i)); 
end 

if b(i) = 0 
mult2(i) =0; 
else 

mult2(i) = (4*d(i)/b(i)); 
end 

co(i) = -0.8572 + 0.5339*mult2(i); % FROM FORMULAE IN REF. 3 

cl(i) = 3.734 - 1.3661*mult2(i); 
c2(i) = -1.7323 + 0.8679*mult2(i); 

Cs(i) = co(i)+(cl(i)*mult(i))+(c2(i)*mult(i) A 2); 
if Cs(i) > 0.918 
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Cs(i) = 0.918 


;% CAPS MAX Cs VALUE AT MAX AREA STA. 


end 

Csdsq(i) = Cs(i)* d(i) A 2; % CALCULATES VECTORS FOR INTEGRATION 

Csdsqx(i) =Cs(i)*d(i) A 2*x(i); 

Csdsqxsq(i)=Cs(i) *d(i) A 2*x(i) A 2; 
end 

% DISTANCE OF EACH ALONG HULL FROM BOW TO STERN 
huldist= [0 9.658 19.315 24.725 30.134 49.91 69.768 89.586 109.403... 

129.22 149.037 168.854 188.671 208.489 228.306 248.123... 

267.94 291.502 315.065 335.628 362.19 385.753 409.315]; 
xi=linspace(0,409,409); 


Csdsqfit=interp 1 (huldist,Csdsq,xi) ;% FIT STATION DATA TO LENGTH OF HULL 

Csdsqxfit=interp 1 (huldist,Csdsqx,xi) ;% FIT STATION DATA TO LENGTH OF HULL 

Csdsqxsqfit=interpl(huldist,Csdsqxsq,xi);% FIT STATION DATA TO LENGTH OF HULL 

int 1 =trapz(Csdsqfit); 

int2=trapz(Csdsqxfit); 

int3=trapz(Csdsqxsqfit); 

m2= (k2*pi)/(L*T A 2)*intl; 

xbar=(int2)/intl; 
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mz=(kprime/k2)*m2; 


kl; 

k2; 

kprime; 

%subplot(21 l),plot(Cs) 
%subplot(212),plot(Csfit) 
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